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Abstract. Microscopic processes on surfaces such as adsorption, desorption, diffusion and re- 
action of interacting particles can be simulated using kinetic Monte Carlo (kMC) algorithms. Even 
though kMC methods are accurate, they are computationally expensive for large-scale systems. 
Hence approximation algorithms are necessary for simulating experimentally observed properties 
and morphologies. One such approximation method stems from the coarse graining of the lattice 
which leads to coarse-grained Monte Carlo (GCMC) methods while Langevin approximations can 
further accelerate the simulations. Moreover, sacrificing fine scale (i.e. microscopic) accuracy, meso- 
scopic deterministic or stochastic partial differential equations (SPDEs) are efficiently applied for 
simulating surface processes. In this paper, we are interested in simulating surface diffusion for pat- 
tern formation applications which is achieved by suitably discretizing the mesoscopic SPDE in space. 
The proposed discretization schemes which are actually Langevin-type approximation models are 
strongly connected with the properties of the underlying interacting particle system. In this direc- 
tion, the key feature of our schemes is that controlled-error estimates are provided at three distinct 
time-scales. Indeed, (a) weak error analysis of mesoscopic observables, (b) asymptotic equivalence of 
action functionals and (c) satisfaction of detailed balance condition, control the error at finite times, 
long times and infinite times, respectively. In this sense, the proposed algorithms provide a "bridge" 
between continuum (S)PDE models and molecular simulations Numerical simulations, which also 
take advantage of acceleration ideas from (S)PDE numerical solutions, validate the theoretical find- 
ings and provide insights to the experimentally observed pattern formation through self-assembly. 
Such phenomena are characterized by a complex energy landscape where the role of noise is critical 
in the emergent behavior of the system. The stochastic fluctuations of the proposed algorithms 
are directly derived from the microscopic model allowing us to explore all experimentally observed 
pattern morphologies starting from a uniform initial state. 

Key words. Interacting particle systems, stochastic (partial) differential equations, Langevin 
approximation, surface diffusion, pattern formation. 

1. Introduction. Surface diffusion of interacting particles as well absorption, 
desorption, reaction, etc. can be accurately simulated using kinetic Monte Carlo 
(kMC) algorithms [1], [2]. In particular, Ising models are set on a lattice and each 
site of the lattice has an order parameter (spin) that describes the presence or not 
of a particle as well its type (Potts models) [3 . Surface diffusion is characterized by 
spin exchange between neighboring sites (Kawasaki dynamics) and depending on the 
rates of the process, kMC evolves the system towards the equilibrium states. However, 
microscopic simulation is computationally expensive when large spatiotemporal scales 
observed in real-life experiments are studied. 

One approach of accelerating the microscopic simulation was developed in a series 
of papers [1], [5], [6] called coarse-grained Monte Carlo (CGMC) method. In CGMC 
setting, the microscopic lattice was coarse-grained and the spins was grouped into 
cells resulted in smaller number of system parameters. Rigorous error analysis was 
performed in [7] and [8] showing that the finite time error is controlled by the interplay 
of the coarsening factor, the temperature and the smoothness of the interaction po- 
tential. Particularly for surface diffusion, it was shown that coarse-graining resulted 
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17 not only in the reduction of the number of system parameters but also in time accel- 

18 eration (square of the coarsening factor faster). In general, CGMC works satisfactory 

19 for long range and mid range interaction lengths, however, it may produce erroneous 

20 results especially for short range interactions, nevertheless, recent variations of the 

21 basic CGMC algorithm have been proposed, trying to overcome this limitation [9]. 

22 Even though CGMC is a powerful tool for accelerating microscopic kMC algorithms, 

23 we are primarily interested in studying pattern formation on surfaces with the ex- 

24 pected patterns having relatively small size. Thus, in order not to lose the necessary 

25 resolution of the patterns, we need to keep the coarsening factor small making CGMC 

26 method a rather inefficient approach. 

27 Another approach to accelerate even further the microscopic simulations is to de- 

28 rive mesoscopic equations by letting the number of interacting particles tend to infin- 

29 ity. Mesoscopic equations for interacting particles are either deterministic or stochas- 

30 tic integro-differential equations. Deterministic PDEs have been used to study nucle- 

31 ation, pattern formation, alloys, etc. [10], [11], [12]. However, thermal fluctuations 

32 (i.e. noise) are important for studying the dynamics thus, more recently, stochastic 

33 PDEs (SPDEs) have been used [13], [H], [H], [H], [17]. For instance, Ostwald ripen- 

34 ing was studied in [17 using an SPDE model. All kinds of numerical schemes such 

35 as finite differences, finite elements, finite volumes as well (pseudo-) spectral methods 

36 have been applied for the discretization of the studied (S)PDEs which leads to a sys- 

37 tem of ODEs in the deterministic case and to a system of SDEs in the stochastic case. 

38 Another technique to derive a system of SDEs that simulates microscopic processes 

39 is by a Langevin approximation of the master equation [18], [19]. Such a Langevin 

40 approximation was derived and studied for surface diffusion and Arrhenius dynamics 

41 in [20^ where it was shown that not only the weak error but also the large deviation 

42 properties of the model are correctly handled. However, in the above studies few or 

43 no care was taken about the exact equilibrium (i.e. invariant) measure of the simulat- 

44 ing process primarily because of the difficulty in satisfying detailed balance condition 

45 (DBC). 

46 In this paper, we develop three different systems of SDEs which serve as approx- 

47 imation models of the microscopic surface diffusion process and additionally satisfy 

48 DBC. The first model is a second-order space-discretization of the mesoscopic SPDE 

49 which is also related with the coarse-grained Langevin (CGL) approximation of [20 . 

50 Even though it is a discretization scheme of the SPDE, we refer to this stochas- 

51 tic model as direct Langevin approximation model (DLM) because its local error is 

52 asymptotically of the same order as the CGL approximation. Furthermore, large 

53 deviation computations show that the action functional between DLM and the micro- 

54 scopic process are asymptotically equivalent thus, rare events and phase transitions 

55 are correctly represented [21 . However, DLM does not satisfy DBC hence its invari- 

56 ant measure so important for determining the equilibrium states or for applications 

57 such as sensitivity analysis of system parameters is not known in general. Nev- 

58 ertheless, the structure of DLM allows the construction of a variant model which 

59 satisfy DBC. Indeed, the second system of SDEs named as perturbed Langevin ap- 

60 proximation model 1 (PLMl) is derived by adding a "correction" term to the drift of 

61 DLM. Then DBC is satisfied and the invariant measure is easily obtained. However, 

62 the "correction" term depends on the coarsening factor hence the cost to be paid is 

63 that the local error is no more as accurate as the local DLM error which results in 

64 perturbed finite time dynamics. 

65 The third model which is named PLM2 eliminates the effect of the "correction" 
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66 term by slightly perturbing the invariant measure in a controlled-error manner. For 

67 Metropolis dynamics, the "correction" term is of diffusion type thus the perturbation 

68 of the invariant measure is an additional higher order term to the entropy. The effect 

69 of this perturbation is that the interacting particle system is simulated at a slightly 

70 different temperature than the original! Similar but more complex perturbation of 

71 the invariant measure is also obtained for the Arrhenius dynamics. Overall, in any of 

72 the proposed models, the error performed either in finite times or in infinite times is 

73 controllable not only to the asymptotic limit but also for any coarsening factor and 

74 actually the interconnection between the finite and infinite time errors as highlighted 

75 by PLMl and PLM2 models is one of the key findings of this paper. 

76 Having derived the stochastic models, the final step in order to simulate and 

77 test them on computers is to discretize the time, too. Since our primal goal is to 

78 highlight the space-discretization properties, we keep the time-discretization as simple 

79 as possible. Thus a simple predictor-corrector (PC) Euler scheme which has 1st order 

80 weak convergence [2T is used. PC Euler which is a two step method can be thought as 

81 a compromise between an explicit and an implicit scheme. Higher order schemes such 

82 as Milstein's or derivative-free Runge-Kutta method could also be applied. However, 

83 higher order schemes are computationally expensive especially for high dimensional 

84 systems such as the studied. 

85 The computational savings of the proposed models compared to the microscopic 

86 system come from many directions. Except for the computational acceleration due 

87 to the coarse-graining which as we already mention is rather limited due to the ap- 

88 plication we are interested in (i.e. pattern formation), there are two other important 

89 acceleration points. The first acceleration is that while in CGMC algorithms only 

90 one particle is allowed to hop between neighboring cells in a time step, in Langevin 

91 approximation more than one "particles" could change their positions on the lattice 

92 in a single time step. The second acceleration stems from the fact that in order to 

93 perform a time step the convolution between the interaction potential and the coarse- 

94 grained lattice configuration is needed to be computed. Convolution can be performed 

95 in Fourier space which results in huge computational savings especially when the in- 

96 teraction potential is long range. This feature is primarily an advantage of spectral 

97 methods which is integrated into our algorithms making eventually the computational 

98 cost independent of the interaction length. 

99 Finally, the proposed Langevin approximation models are applied to the study of 

100 pattern formation phenomena. Such phenomena are characterized by a complex en- 

101 ergy landscape [24] , [25] , [26] where the role of noise is critical in the emergent behav- 

102 ior of the system. The stochastic fluctuations of the proposed algorithms are directly 

103 derived from the microscopic model, they allow us to systematically explore all exper- 

104 imentally observed pattern morphologies through a self-assembly mechanism, starting 

105 from a uniform initial state (non-equilibrium dynamics). Indeed, using Morse- type 

106 interaction potential, which is an attractive/repulsive potential, at various parameter 

107 regimes, we were able to reproduce the experimentally observed 2D images shown 

108 in [2T. Moreover, we study different versions of Morse potential so as to reveal the 

109 importance of stochastic fluctuations not incorporated in other analysis tools such 
no as linear stability analysis or deterministic PDEs which are usually trapped in local 

111 minima of the complex energy landscape. 

112 The organization of the paper is as follows. Section [2] discusses the microscopic 

113 Ising formulation for surface diffusion as well the coarse-grained model for Metropolis 

114 and Arrhenius dynamics. Langevin approximation and mesoscopic SPDEs for both 
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115 dynamics are also presented in the same Section. In Section [3) the proposed SDE 

116 models are presented and their approximation properties are derived while in Section [4] 

117 pattern formation phenomena are observed and studied. Finally, Section [5] concludes 

118 the paper and suggests further directions of future work. 

119 2. Background. Let us begin with the presentation of the microscopic model 

120 and continue with its coarse-grained analog. Two different dynamics namely Metropo- 

121 lis and Arrhenius are considered. Then the derivation of CGL approximation model is 

122 reviewed and finally the mesoscopic SPDEs which one way to be obtained is through 

123 taking the limit of the coarse-grained model [4 are given. Fig. |2.1[ a) schematically 

124 depicts the position in space and time scales of the revised models while Fig. |2.1[ b) 

125 shows the actual lattices of various models discussed in the following Sections. Please 

126 note that our interest in this paper lies both in microscopic and in mesoscopic scales. 

2.1. Microsco pic Model. Consider a finite, periodic, d-dimensional, fine lattice 



128 (left drawing in Fig. 2.1 ^b)) defined by Cn '= j^^^ HP^ where is the size of the 

129 lattice site while is the total number of sites of the lattice. At each lattice site 

130 X G Cn^ an order parameter -usually referred as spin- is allowed to take two values 

131 describing vacant and 1 describing occupied. On the fine lattice a spin configuration 

132 is defined as a := {(j{x) G {0, 1} : x G Cn} and it is an element of the configuration 

133 space S := {0,1}^^. 

134 The energy of the system evaluated at a is given by the Hamiltonian 

H{a):=-\ J{x-y)a{x)a{y)+ h{x)a{x) (2.1) 

x,y e Cn x^Cn 
y ^ X 

135 where J(-) is the interaction potential between the sites while h{-) is the external field 

136 applied to the system. Note that the interaction potential has radial symmetry and it 

137 is appropriately scaled so as the derived mesoscopic limit is well-defined. Moreover, 

138 interaction potential has support in [— ^, ^]^, thus, its interaction length is L sites. 

139 Equilibrium states (i.e. invariant measure) of the model at inverse temperature f3 is 

140 described by the Gibbs measure given by 

^NAd':^) = -^e-^"<^^^Pida) (2.2) 

141 where ^Ar,/3 is the normalization factor that makes /iAr,/3 a measure while P{da) is the 

142 prior measure defined as a product of independent Bernoulli random variables one for 

143 each lattice site. 

144 Surface diffusion is simulated as spontaneous spin exchange between two neigh- 

145 boring sites, x, y with the restriction that every site cannot contain more than one 

146 particle (exclusion principle). Two different spin exchange dynamics which satisfy 

147 the detailed balance condition are considered. The first surface diffusion dynamics 

148 is Metropolis and its exchange rate is defined for two by neighboring sites, x, y at 

149 configuration a as ^ 

c{x,y,a) := c^o exp(/3 min{0, {a{x) - a{y)){U{x,a) - U{y,a))}) (2.3) 



150 where do is the diffusion rate which depends on physical properties of the surface 

151 while U{x^a) is the potential of the site x given that the current configuration is a 
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Hierarchical (Multiscale) Modeling 
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(a) Due to wide range of characteristic length scales and characteris- 
tic time scales, several models and simulation algorithms have been 
developed in the literature. The proposed Langevin models provide a 
"bridge" between the continuum models and the microscopic processes. 



Microscopic Coarse-grained Langevin 

Lattice l-attice l-attice 




(b) Various lattices at different scales in space. Notice that both coarse- 
grained and Langevin lattices have the same and known spatial scale 
while their order parameters take discrete and continuous values, re- 
spectively. 

Fig. 2.1: Hierarchical modeling at different time and space scales. 



152 defined as 

U{x,a):= J2 J{x-y)a{y)-h{x) (2.4) 

y G Cn 

153 Despite using Metropolis dynamics in many studies, a more natural and possibly 

154 more appropriate description of the finite time surface diffusion dynamics is Arrhenius 

155 dynamics. In Arrhenius dynamics, a spin exchange is performed when the activation 

156 energy is above an energy barrier which depends on the properties of the potential en- 

157 ergy of the surface [28] , [29] . Arrhenius dynamics for spin exchange (surface diffusion) 
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158 between two neighboring sites, x^y is given by 



c{x, y, a) := do{l - a(x))a(^)e-^(^°+^(^'")) + doa{x){l - a(y))e-^^^°+^^^'")) (2.5) 

159 where do and Uq are the diffusion rate and the energy barrier of the surface, respec- 

160 tively, and depend on the physical properties of the diffusion process while U{x^g) is 

161 as before the potential of the site x. Thus, a continuous-time jump Markov process 

162 {crt}t>o on L'^(i;;R) is defined with generator 



174 
175 



182 
183 



E[/K)H=/:/(<t)= c{x,y,a)[f{a^'''y^)-f{a)) (2.6) 



d 
~dt 

X ^ y 



163 for any test function / G L^(i;;R). Please note that a*^^'^^ denotes the new configu- 

164 ration of the lattice after one spin exchange between neighboring sites x and y while 

165 test function / also called observable is typically independent of the size of the lattice. 

166 A special class of observables called mesoscopic plays a crucial role in the proofs of 

167 the approximation theorems in [7], [30]. 

168 2.2. Coarse-Grained (GC) Model. The coarse-graining of the microscopic 
system is performed by grouping the sites of the microscopic lattice into cells. Each 
cell is denoted as Ck with size \Ck\ = where q is the coarsening factor at each 
dimension while k G Cm where Cm •= HP^ lattice (middle drawing 



in Fig. 2.1 ^b)). Obviously, the size of the CG lattice is with m = N/q. On the 



CG lattice, a CG variable is defined for the kih. cell by 

r]t(k) := ^ (Jt{x), keCm (2.7) 

xeCk 

thus a new continuous-time jump Markov process {r]t}t>o is defined. In what follows, 
our primal interest is concentrated on the averaged coarse-grained variables defined 

176 as 

Uk) := keC^ (2.8) 

177 which are elements of the configuration space Hg,rn = {0^ •••^ 1}^"^. 

178 As in the microscopic formulation, averaged CG process has Hamiltonian, poten- 

179 tial, rate (dynamics) and invariant measure which are approximations of the respective 

180 microscopic quantities. The Hamiltonian of the averaged CG process is given by 

H{fi):=-^ E Jik-l)fikm+ T.CHk) + ^)fik (2.9) 



where J(-) is the coarse-grained interaction potential given by 

j{k-i):={ , r ; (2.10) 



hY.lfcT-J^^-y) for k^l 

Eltlco'^i^-y) for k = l 



Equilibrium states of the averaged CG variables at inverse temperature /? has invariant 
measure given by 

l^,,mAdv) = ir^e-'l"^"^'^^P,,ra{df]) (2.11) 
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184 where Zq^rn,(3 is the normahzation factor that makes iJi{dff)q^ra,p a measure while 

185 Pq,m{dv) is the prior measure defined as a product of binomial random variables 

186 one for each coarse cell. 

187 The rate of the averaged CG process to jump a particle from a cell /c to a neigh- 

188 boring cell denoted by Ck^i{ff)^ is given for Metropolis dynamics by [5] 

CkA^) '•= doq%{l - m)exp {/3 mm{0,U{l,f]) - U(k,f])}) (2.12) 

189 where 

u{k, fj) := ^(^ - 0^(0 - cm + m) (2.13) 

190 is the CG potential of the kth cell. On the other hand, the exchange rate of a particle 

191 between two neighboring cells /c, / is given for Arrhenius dynamics by [5j 

CkAv) ■■= doq%{l - r?Oe-''(''°+^(''*^» (2.14) 

192 Thus, the generator of the averaged CG variables, {f]t}t>o^ is 

^E[f{f),M = cf{f])= ^ ckAmiv+^im-Skm-fm (2.15) 

193 for any test function / G R). 

194 Finally, the weak error analysis between the microscopic process and the CG 

195 process performed in [7] uses consistency with the backward Kolmogorov equation 

dtW^jCw = , t<T ^ ^ 

196 which corresponds to the master equation for expected values w{z^t) = 'E[f{f]T)\f]t = 

197 z]. Thus, using observables with bounded derivatives and Kolmogorov consistency, it 

198 was rigorously shown that the weak error between the microscopic process and the 

199 CG process is of order 0((^)^) which is affordable for mid and long range interaction 

200 potentials {L » 1). 

201 Remark: The computational acceleration of the CGMC algorithm for simulating 

202 surface diffusion processes stems not only from the reduced number of parameters by 

203 a factor of but also from the time acceleration by a factor of q'^^ [6 . Intuitively, the 

204 time-acceleration can be understood by the fact that one event in the CG simulation 

205 is the jump of a particle from a cell to a neighborhood cell while in microscopic 

206 simulation the same event can be a (possibly long) sequence of jumps. 
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2.3. Coarse-Grained Langevin (CGL) Approximation. Generally in Langevin 
methods, the microscopic process is approximated by a process driven by a system 
of SDEs [18 , [31 . For surface diffusion particularly, Langevin approximation for 
the coarse-grained model was recently derived in [20 . Concentrating for notational 
simplicity in ID, the CG Langevin SDE system is given by 

dVk = ak{f})dt + hkMdWu k e Cm (2.17) 

where fj = {fjk ' k G Crn} is the S DE process set on the configuration space Hq^rn = 
[0, 1]^"^ (see right drawing in Fig. 2.1 [b)) while a{fj) = {ak{fj) : k G Cm} and b{fj) = 
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214 {bk,i{Tj) ' k^l e Cm} are the drift vector and the diffusion matrix of the SDE process, 

215 respectively. The generator of this process is defined for an arbitrary test function 

216 /gL^(H,,^;R) as 



217 
218 
219 



In order to estimate the drift and diffusion terms, the weak global error between 
the CG process and CGL process is minimized. Thus, defining for a mesoscopic 
observabl^ /, the expected value w{z^t) = ¥.[f {f}T)\f}t = weak error is written as 

E[f{fjT)] - nfiVr)] = m[f{VT)\VT = Vt]] - m[f{VT)\fjo = Vo]] = 

E[w{vT. T)] - E[w{vo. 0)] = ^ E[Cw{v) + dtw{f])]dt = ^^.ig) 

E[Cw{f]) - Cw{f])]dt = [ E[eioc{w)]dt 



220 where the third equation is the martingale property w hile th e fourth one uses the 

221 backward equation for jC [30 . Moreover, according to (2.19), the local error for a 

222 mesoscopic observable /, eiodf)^ can be defined on the difference of the generators of 

223 the two processes as 

eiocif) = Cf{n) - Cf{fj) 

= E ck,imm+-m)-Skm-w{fj)) 

k,iec^ ^ (2.20) 

224 2.3.1. Weak Error Analysis. By applying Taylor series expansion for f{f] + 

225 ^{Si{k) — Sk{l))) and appropriately choosing the drift and diffusion terms so as to 

226 eliminate the first and second order of the expansion, it was obtained in [20 that the 

227 kth element of the drift vector is 

akiv) = ^[c/c+i,/c(^) - Cfc,/c+i(^) + Ck-i,k{v) - Ck,k-i{v)] (2.21) 

228 while the non-zero elements of the diffusion matrix are 

h,k{v) = ^\Jck+i,k(f})^Ck^k+i(r}) ^2 22^ 

h+i,k{ri) = -h,k{rj) 

229 Thus the formal local error between CG process and CGL approximation process is 

eiociw) = 0(4) X 0{ck,i) = 0{\) (2.23) 



^ A mesoscopic observable is a function whose derivatives -in this particular case up to third 
order \30\- are bounded and the bounds are independent of the dimension (i.e. the size of the CG 
lattice). 



Therefore, based on the above approximation, finite time global weak error be- 
tween CG process and its CGL approximation could be rigorously obtained for meso- 
scopic observables by using again Kolmogorov consistency of the backward equation 
and Bernstein-type bound estimates for the derivatives of w{z,t). Indeed, it was 
shown in [30 that the weak error for mesoscopic observables is 

E[/(f?T)]-E[/(f?T)]=0(^) (2.24) 

when absorption/desorption processes were considered and we expect the same result 
is true for diffusion processes. 

2.4. Mesoscopic SPDE Limit and LDP. In this Section, we review meso- 
scopic evolution equations arising in surface processes derived from the microscopic 
stochastic models presented above. In general there are two families of mesoscopic 
equations depending on the presence of stochasticity. Here we concentrate on the 
stochastic integro-differential equations for Metropolis and Arrhenius dynamics. Both 
dynamics can be written as a constrained gradient fiow equation plus a multiplica- 
tive stochastic term with different mobilities. Indeed, the unified stochastic mass- 
conserved equation (SPDE) is given formally by [11], [15] 



dtp 



245 where t) is the zero lattice-size limit of the empirical measure of the particles 

246 which evolves slowly similar to a density while E[-] is the Lyapunov functional (free 

247 energy functional) of the deterministic mesoscopic equation given by 

E[p] = -^ J J J{x - x')p{x)p{x')dxdx' ^ 13 j h{x)p{x)dx^ J R{p{x))dx (2.26) 

248 where J(-) and h{-) are continuous versions of the interaction potential and external 

249 field, respectively, while R{-) is the entropy of the system given by 

R{p) = p\og{p) + {l-p) log(l - p) (2.27) 

250 L[p] is the mobility of the equation which determines the dynamics of the system 

251 while W{x^t) is space-time white noise. The invariant measure for the equilibrium 



252 states of the solution of (2.25) is given formally by |2j 



Mdp) = -^e-'^'^^PUp (2.28) 
253 A formal approach to derive the above invariant measure is to take the zero lattice- 



254 size limit of the CG invariant measure given by (2.11). Indeed, another way to write 



255 down (2.11) is to expand the binomial prior distribution using Sterling's formula [4 . 



256 Then the invariant measure is written as 



where 



,,^^Adv) = _^e-^-^(^(^)+i^'5(^^+°(^^» (2.29) 



E{fj) = \m{v)+m] (2.30) 



258 is a disc rete version of the Lyapunov functional while H{-) and R{-) are the Hamilto- 
nian in i2M and the discrete entropy of the system (i.e. R{f]) = ^keCm^^^ ^^^iVk) 



259 



269 



270 



276 
277 
278 



280 



260 {l — r]k) log(l — r]k)])^ respectively. The additional term in (2.29) is the primal remain- 

261 der of the Sterling's expansion which equals to G{f]) = ^TYlkeCm — Vk))- 

262 Notice that the additional term, G(-), may be significant when coarsening factor, 

263 takes small values, however, in the zero lattice-size limit the only term that survives 

264 is the Lyapunov functional, E{-). 

265 The mobility for Metropolis dynamics equals to L[p\ = dop{l — p)^ thus, the SPDE 

266 for Metropolis dynamics is given by 

dtP = V- {doiVp - /3p{l - p)V{J *p))} + -^V ■ { V^dopil - p)w} (2.31) 

267 where * denotes convolution. For Arrhenius dynamics, the mobility is more complex 

268 and it is given by L[p\ = dop{l — p) exp(— /3(/7o + J * p)), thus, the mesoscopic SPDE 
for this case is 

dtP = V- {d^ exp(-/3J * p){Vp - pp(l - p)V{J * p))} 

+ V • I ^2d^p(l-p)eM-PJ^PW^ ^^'^^^ 



where d^ = d^e 



-/3t/o 



271 Finally, SPDEs such as (2.25) are generally ill-behaved mathematical objects 

272 especially in high dimensions and they are usually treated in a formal way as here. 

273 Nevertheless, an indirect yet rigorous analysis could be carried out for SPDEs using 



274 the theory of Large Deviations (LD) ^32^. Indeed, SPDE (2.25) is related with the 



275 action functional for the microscopic process obtained by taking the hydrodynamic 
limit. For exchange dynamics with exclusion principle and long range interaction 
potential, it was shown in [33 that the action functional for an absolutely continuous 
function ^ : [0, 1]^ x [0,T] ^ R equals to 



^ot(^)=/ / L[^{VHfdxdt (2.33) 
Jo io 



where H solves 



= V • [m{^^^^ - /?V(J * * + + 2V • {mWH} (2.34) 

which is the second order backward PDF of ( |2.25[ ). Intuitively, the action functional 

281 Sqt{^) assigns a probability to the event p that follows the path ^ which can be 

282 formally stated by the following asymptotic formula 

P{v{p, ^)<5}r^ ^-n-^SotW (2.35) 

283 for suitably chosen ^, e > where u is a metric in a proper function space that mea- 

284 sures the distance between p and ^. Further details on LD theory can be found in 

285 Section 13.41 

286 Remark: Even though, mesoscopic models -either deterministic or stochastic- are 

287 computationally tractable compared to microscopic or even CG models they lack of 

288 some interesting properties. For instance, due to the limiting process, the actual 

289 length-scale of the system is not obvious. Moreover, the space discretization is not a 

290 trivial issue especially for the stochastic case since the properties of the discrete and 

291 the continuous models may be totally different. These facts will be highlighted in the 

292 following Sections. 
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3. Langevin-type Approximation Models. As it was reviewed in Section 2.3 
classical Langevin models are derived as approximations of the atomistic processes by 
formally minimizing the local error between the microscopic process and the SDE pro- 
cess. In connection with Fig. |2.1[ a), Langevin approximation is an approach which 
translates the atomistic processes from the microscopic level to the coarser mesoscopic 
level. In this Section, we proceed in the opposite direction (i.e. from mesoscopic to 
microscopic level) and derive three Langevin-type models from mesoscopic equations 
for the simulation of surface diffusion processes which additionally to the properties 
of the classical Langevin approximation they satisfy - actually two of them- detailed 
balance condition (DBC). Eventually our goal is to control the error of the derived 
approximations at three different time-scales which are 

a. Finite times through weak error estimates between the microscopic process 
and the derived models. 

b. Long times and phase transitions through LD theory and asymptotic equiv- 
alence of the rate (action) functionals. 

c. Infinite times through the knowledge of the invariant measure of the derived 
approximation process. 

To begin, the first proposed model is a 2nd order space discretization of the 
mesoscopic SPDE. We refer to it as Direct Langevin approximation model (DLM) 
because the local error between DLM and CGL of [20 is of order (-^ ) for the drift 



term while it is of order O(^) for the diffusion term (see Section 3.1.1) which are 
considered negligible. Yet, as in CGL approximation, the DBC is not satisfied for 
DLM thus the invariant measure of the stochastic process is not known. By adding 
a "correction" term to the drift, the second model referred to as perturbed Langevin 
approximation model 1 (PLMl) is defined. For this variant, DBC is satisfied and 
the invariant measure is a discrete version of the continuous invariant measure given 
by (2.28). However, the finite time dynamics of PLMl are perturbed due to the 
additional "correction" term. The third model referred to as PLM2 tries to overcome 
the induced error at the dynamics by adding a perturbation term to the invariant 
measure. An appropriate choice of the perturbation term leads to the elimination of 
the "correction" term of PLMl restoring the accuracy of the finite time dynamics. 



Table |3TT] summarizes the properties of CGL approximation as well the properties of 
the three proposed models which we will derive in the remaining of this Section. 





Weak Error of order O(^) 


LD Theory 


Invariant Measure 


CGL 


Yes 


Yes 


No 


DLM 


Yes 


Yes 


No 


PLMl 


No 


Yes 


Yes 


PLM2 


Yes 


Yes 


Yes 



Table 3.1: Summary of the properties of the derived diffusion models at different time 
scales. Note that for adsorption/desorption processes the answer to the LD Theory 
column is 'No' [34]. [20]. 



326 Before starting presenting the proposed models, we make the following simplifi- 

327 cations. Without loss of generality we concentrate on the ID case. We revisit the 

328 general d-dimensional case in Section [4] where the details of the numerical implemen- 

329 tation are given. Moreover, external field is assumed to be zero without this being a 

330 restriction to the final outcome. 



11 



332 
333 
334 
335 



344 
345 



350 
351 
352 



358 
359 
360 



3.1. Direct Langevin Approximation Model (DLM). The first approxima- 
tion model is a straigiit forward second-order, finite-difference, mass-conserved space- 
discretization of the mesoscopic SPDE. The discretized density vector is denoted by 
p = {pk : k e Cm}' Then for the kth density element, a stochastic differential 
equation is defined by 



dpk = Uk{p)dt + ^ Vk^i[p)dWi, k e Cr, 



where 
Uk{p) = 



dE{p) dE{p) 
dpk+i dpk 



dt 



{Lk{p) + Lk-i{p)) 



(3.1) 



dE{p) dE{p) - 

dpk dpk-i 
(3.2 



is the kth element of the drift vector. Note that E{p) is the discrete free energy 
functional given by (2.30) while L}^{p) is the discrete version of the mobility. For 
Metropolis dynamics, the mobility is given by 



Lk{p) = dopk{l - pk) 



(3.3) 



which depends only on the /cth density parameter pk while the mobility for Arrhenius 
dynamics is given by 

Lk{p)=d^Pk{l-pk)e-^^^^^'^ , (3.4) 

which depends not only on pk but also on the neighboring density variables through 
the potential U{k^p). The non-zero elements of the diffusion matrix are 



^k,k{p) 
^fe+l,fe(p) 



(3.5) 



-(Lfc+l(p) +L/e(p)) 

-Vk,k{p) . 

Hence the covariance matrix (i.e. square matrix of the diffusion matrix) is a tridiagonal 
matrix with non-zero elements 

1, 



{vv^)k,k{p) = ^[Lk+i{p) 
{vv^)k±i,k{p) = 



Lk-i{p)^2Lk{p)] 
[Lk±i{p) ^ Lk{p)] . 



(3.6) 



It is noteworthy that the scaling of the noise in (3.5) is ^ which is different 
from the scaling -^i^ of the mesoscopic SPDE (2.25). The reason is that in order to 
relate the process generated from (3.1) with the CG process or the CGL process (i.e 
Pk ^ fjk ^ fjk) the appropriate scaling for the stochastic term is ^ as the following 
subsection reveals. Linked with Fig. |2.1[ a), different scalings of the noise result in 
models with different positions at the mesoscopic level. Typically, when zooming into 
the atomistic details is performed, the power of the noise is increased while when 
zoom out is performed the noise is faded out. 

Additionally, the existence of a Lyapunov functional is usually crucial for the 
study of an (S)PDE either theoretically or numerically. In (3.1), if the noise is can- 
celled out then E{p) is a discrete Lyapunov functional since it is decreasing over time 
(see Appendix [b|. Of course, when noise is present Lyapunov functional may increase 
due to stochastic fluctuations nevertheless on average it decreases. Next we proceed 
with the properties that relates the process driven by (3.1) with the CG and CGL 
processes. 
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3.1.1. Weak Error Analysis. The estimation of the finite-time weak error be- 
tween the DLM process, p^, and the CG process, fjt^ uses as an auxihary intermediate 
step the CGL process, f]t. Indeed, the weak error for a suitable mesoscopic observable, 
/, can be written as 

E[/(^?t)] - E[/(pt)] = nfiVT)] - nfiVT)] + nfivr)] - E[/(pt)] (3.7) 
and at least formally it was shown in [2F and briefiy reviewed in Section |2.3| that 



HfiVr)] - IE[/(7)t)] =0(4). On the other hand, the local error between the CGL 



process and the DLM process defined in (3.1) is given by 



Cf{p)-Mf{p)= Y,[au{p)-uu{p)]^^-\ [ibb^)M{p)-{vv^)ki{p)]j^^ 



dpkdpi 
(3.8) 



where Ai is the generator of the process driven by (3.1) given by 

' dpkdpi 



keCr, 



(3.9) 



369 
370 
371 



for any test function / G L^(Hg^^;R). 

It is straightforward to compute (see Appendix [A| that the drift term has the 
following formal asymptotic expansion 



^fe(p) = { Lk{p) 



dxp{xk) 



p{xk){l - p{xk)) 



+ 0(-^) 



(3.10) 



372 
373 
374 
375 



where p{xk) = p{xk^t) is the continuous space density function at position = 
— , /c = 0, ...,m — 1 and it should not be confused with the DLM process, p^, which 
is discrete in space. Similarly, the weak asymptotic formula for the covariance matrix 
of the diffusion for two test functions ^i(x) and 02 (^) is given by 



^Ik {P)^i i^j ) . Yl (P) ^2 {Xi ) 



2 

qm 



i,i 



dt 



[ L[p{x)]dx(l)i{x)dx(l)2{x)dx^0{^) 
J 



(3.11) 



376 
377 
378 



The same asymptotic expressions have been derived for CGL approximation in [20 . 
Moreover, applying the time rescaling t m?t suggested by the above asymptotics 
to both DLM and CGL processes, it is allowed to formally write that 



Uk{r}) = o.k{r}) + O(^) 



(3.12) 



379 
380 
381 



where a(-) is the drift vector of the CGL process given by (2.21). Similarly, having 
in mind that Brownian motion scales as Wrn2t = m^^' straightforward to show 
that 



{vv'^)k,k{v) = {bb^)kM + 0{ — ) 

qm 

{vv^)k±iMv) = {bb^)k±iAv)^0{^) 

qm 
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(3.13) 



382 where b{-) is the diffusion matrix of the CGL process. Thus substituting (3.12) and 

383 (3.13) into (3.8) we derive at least formahy that 



393 
394 



395 
396 
397 
398 
399 
400 

401 
402 
403 
404 



405 
406 
407 
408 
409 



Cfip)-Mf(p) = Oi—) 
qm 



(3.14) 



and using the same arguments presented in and briefly reviewed in Section 2.3 
the weak error could be rigorously proved to have the same O(^) order of error 
Finally, notice that q « m hence the weak error between the CG process and the 
DLM process is of order O(^). 

3.1.2. Is DBC satisfied? A guess for the invariant measure of the DLM process 
could be 



1 



(3.15) 



390 which is a discrete version of (2.28). However, this guess is not correct because the 



391 operator JVl (i.e. the generator) is not self-adjoint ^ with respect to the 

392 measure \i. Indeed, we compute (see Appendix [B|) that 



2(7 J ^ 



Ck{p) 



dg r df 
opk opk 



p{dp) (3.16) 



where < •, • >L^fi) 
measure p while 



denotes the inner product between two functions with respect to 



Ckip) = 



dL 



k+l 



dL 



k-l 



dpk dpk dpk dpk+i dpk+i dpk-i dpk-i 



(3.17) 



is an interference term which depends only on the mobility of the process. 
Remark: For the case where the mobility is constant (additive noise) or even linear 
then Ck{p) = for all k thus DBC is satisfied and p{dp) is the invariant measure of 
the process. However, the mobility of both Metropolis and Arrhenius dynamics which 
partially reflects the exclusion principle of the microscopic process are more complex 
hence the invariant measure is not known explicitly. 

3.2. Perturbed Langevin Model 1: Satisfying the DBC. The second ap- 
proximation model (PLMl) is obtained by adding a "correction" term to the drift 
which cancels the interference term in (3.16). Thus, the kth element of the density 
vector of PLMl is given by 

dpk = (^Uk{p) + ^Ckip)^ dt^ ^kAp)dWi, keCm (3.18) 

which is obtained from DML with a perturbation of order O(^) to the drift. 

Proposition 3.1. The stochastic process driven by (S.IS^ satisfies the DBC and 
its invariant measure is p{dp) given in 1^3.15 ). 

Proof. The generator of the new process denoted by M. is written for a test 
function / as 



Mfip) = Mfip) + ^^ J2 Ckip) 



keCr, 



dl 
'dpk 



(3.19) 
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410 hence using (3.16 ) which has been derived in Appendix[B]it is straightforward to show 

411 that 



< Mf,g >m^)=< f.Mg >m^,) (3.20) 

412 □ 

413 3.2.1. Weak Error Analysis. Due to the "correction" term added to the drift, 

414 the finite time dynamics of PLMl are perturbed. Indeed, the local error between the 

415 PLMl process and the DLM process for a test function /, which is defined as the 

416 difference of the two processes' generators (see Section 2.3), is 

Mfip) - Mfip) = Ck-^^ = 0( J) (3.21) 

417 Hence, the weak error between the CG process and the PLMl process is expected 

418 to be of order O(^) which is worse than the weak error between the CG process 

419 and the DML process. Overall, the cost paid for constructing a model with known 

420 invariant measure is to introduce error at finite times. Thus, in order to gain better 

421 understanding of the induced error, lets compute explicitly as well asymptotically the 

422 added "correction" term. 

423 3.2.2. "Correction" Term Asymptotics. For Metropolis dynamics, the "cor- 

424 rection" term is twice the discrete Laplacian of the density thus its asymptotic is given 

425 by 

Ckip) = 2[pfe+i + pk-i - 2pk] = ^d^^p{xk) + 0{\) (3.22) 

426 Interestingly, the Laplacian of the density is also obtained asymptotically from the 

427 entropy term of the free energy functional (see (2.31)). Similarly, the "correction" 

428 term for the more complex Arrhenius dynamics is given by 

\pk[l-pk) J 

^ - /?(J(0) + JW)) Lu+^{p) - ( ^-^P'^-' - /5(J-(0) + J{1))] Lu-,{p) 

= { ( p(.l)(wl)) - + ^"(^») ^^('')} + ^(i) 

(3.23) 

429 where the last equation is its asymptotic expansion. However another less accurate 

430 yet more manageable asymptotic expansion for the Arrhenius "correction" term is 

431 needed which is given by (see Appendix [A|) 



(3.24) 



p{xk){l - p{xk)) J q^m^^ 

where 7 = (Xlz^o 1 ^(0) ^ constant. 
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433 3.3. Perturbed Langevin Model 2: Perturbing the invariant measure. 

434 Previous subsection motivates us to suggest a second variant of DLM with perturbed 

435 invariant measure which is able to ehminate the "correction" term from the drift. 

436 Hence, the price to be paid for correcting the finite time dynamics is a controhed- 

437 error approximation of the invariant measure. To proceed, the third approximation 

438 model (PLM2) is derived by assuming that the invariant measure of the DLM process 

439 is a perturbed version of ii{dp). Indeed, assuming that the (perturbed) invariant 

440 measure is given by 



446 
447 
448 
449 
450 

451 
452 
453 
454 
455 
456 
457 
458 
459 



jl{dp) 



1 



W dpk 



(3.25) 



441 where P(-) is a function to be specified, then, the following computation similar to 

442 (3.16) is obtained for the generator M. of DLM 



dg_ 
dpk 



f 



443 where Ck{p) is given in (3.17) while 

dP 



Pk{p) = (Lfc+i(p) + Lfc(p)) 
2 



_dpk+i 



dP 

dpk 



{Lk{p) + Lk-i{p)) 



df 

^ — ^ 
opk 



dP 



dP 



dpk dpk- 



tOx < ft 



jl{dp) 
(3.26) 

(3.27) 



444 is the interference term due to the perturbation of the invariant measure. Then 

445 PLM2 is defined for the /cth density variable by 



dpk = (^k{p) + dt^ ^ Vk^i{p)dWu 



(3.28) 



where Ck{p) = Ck{p) — Pk{p) is the new "correction" term. Similarly, to the previous 
model, PLM2 was an explicitly known invariant measure. 

Proposition 3.2. The stochastic process driven by (3.2^ satisfies the DEC and 
its invariant measure is p{dp) given in (3.25). 

The proof is omitted because it is similar to the proof for PLMl. 

3.3.1. Weak Error Analysis. Choos ing a ppropriately the perturbation term, 
it is possible to make Ck{p) negligible, e.g. (3.30). Eliminating the "correction" term 
implies that the drift term is not anymore perturbed and the finite time dynamics are 
again as accurate as the DLM dynamics. The choice of the appropriate perturbation 
of the invariant measure is inspired by the asymptotic expansions of the interference 



term (3.17) and the invariant measure perturbation (3.27). As already stated, the 



asymptotic of the interference term for Metropolis dynamics is the Laplacian of the 
density hence a suitable choice for the perturbation term is the entropy of the system. 
Indeed, if we set 



^(P) = ^^^(Pk) + (1 - Pk) log(l - pk)] 



(3.29) 



460 then it is obtained asymptotically that Ck{p) = 0{^). Hence the local error between 

461 the time rescaled DML process and the time rescaled PLM2 process for any test 
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462 function, /, is 



Mf{p)-Mf{p) = 0{^) (3.30) 

463 Interestingly, the perturbation term, of the invariant measure for Metropohs 

464 dynamics is the entropy of the system. This imphes an increase of the temperature 

465 of the system at equihbrium from /3 to /3{1 + ^)! Moreover, PLM2 can be thought 

466 as a space discretization of the SPDE ( |2.25| ) since it differs from DLM, which is a 

467 straightforward space discretization of the same SPDE, by a term which has order 

468 less than the order of the discretization. Consequently, it can be stated that numerical 

469 simulations of the discretized process -possibly any discretized process- are performed 

470 at a different (of order 0(^)) temperature than they were initially designed. 

471 For Arrhenius dynamics, the derivation of the perturbed term is more difficult 



472 since the asymptotic expansion given by (3.24) is more complicated. Nevertheless, if 

473 we set 

Pip) = E + (1 - - ^^)] - + j{i))H{p) 

kec^ (3.31) 

~ ^ E log(P/c) - (1 - Pk) log(l - Pk)] 

474 then the asymptotic order of the "correction" term becomes Ck{p) = 0( ^2^4 )• Hence 

475 the local error between the DLM process and the PLM2 process for Arrhenius dy- 

476 namics is given by 

Mf{p) - Mf{p) = O(^) (3.32) 

477 where L is the interaction potential length. Finally, notice that for Arrhenius dynam- 

478 ics both Hamiltonian and entropy terms are perturbed and there is no straightforward 

479 physical interpretation of the perturbation as there was for the Metropolis case. 



Remark: Comparing the perturbations terms (3.29) for Metropolis dynamics and 



481 (3.31) for Arrh enius dynamics with the additional term G(-) in the invariant measure 



482 of CG process (2.29), we observe that they have the same order, O(^), but the actual 



483 functions are different. Of course, this is not a surprise since the former depends on 

484 the mobility (i.e. dynamics) while the latter depends on the prior distribution of the 

485 process. 

486 3.4. Large Deviation and Action FunctionaL It was shown firstly by Hanggi 

487 et al. [34] that Langevin approximation may have different behavioi]^ at long times 

488 compared to the microscopic process. This is established by showing the asymptotic 

489 non-equivalence of the large deviations of the derived models and the microscopic 

490 process as defined by their action functionals. Hence, apart from the local error, 

491 we are interested in the long time behavior of the derived approximation processes 

492 including rare events and phase transitions. It was shown in [21 , where an action 

493 functional for the mean field Ising model was derived, that the asymptotic equivalence 

494 of the action functionals between two processes implies that the processes have similar 



^See the remark at the end of this subsection for such an example. 
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dynamical properties and particularly they have the same probability of rare events 
and exit times. 

In this Section, a time dependent action functional is derived for the DLM. Similar 
computations for the variants of DLM give the same asymptotic behavior thus they 
are omitted. We show that the action functional for DLM is asymptotically equivalent 
to the action functional derived in [33 and briefly revised at the end of the Section 2.4 



where large deviations for a system of long range interactions that models diffusion 
of interacting particles was studied. The results in [33] where an extension of the 
large deviation results in [35 in which Kawasaki dynamics (i.e. diffusion) for short 
range interactions was examined. Since DLM is a space discretization of the SPDE 
(a.k.a. the action functional of the microscopic model) it is straightforward to show 
the asymptotic equivalence of the action functionals. Nevertheless, we present the 
detailed derivation for completeness. 

In order to recover the action functional we have to identify a small parameter 
which will be sent to zero. In our case, the small parameter is the spacing of the 
discretization, ^, or, in the context of coarse graining the size of a cell. Then for any 
absolutely continuous functions ^ : [0,1] x [0,T] R and G : [0, 1] ^ R the rate 
function is given by 

^S^(^)= / A'^(^,^t)^^t (3.33) 

^0 



where 



A^(^,^,) = sup|< g.dt"^ -m^u{ilj) >p < g, qm^vv^ {ilj)g >i2 } , (3.34) 



514 while g = {gk = G{xk)} e M^, similarly V^(t) = {t/jkit) = ^(t, Xk) e and < •, • >^2 

515 is the usual P inner product. 
Using the asymptotic approximations ( 3.1Q| ) and ( 3.11[ ) (i.e. the drift and the 



diffusion of the SPDE) it is straightforward to show that as m ^ oo 

<g,^t- mM4'),9 > =< G, - - /3d,{J * *))| >,2 

(3.35) 

and 

< g, qm^{^)g >,2^< d^G, L[*]a^G (3.36) 
Thus as m — > oo the asymptotic hmit for A"*(5', ^t) is 

A(^, ^t) = sup Gd, - f^9M * ^))} - 1^ mid^Gfdx 

(3.37) 

Using F-convergence arguments and the arguments in [33j, a rigorous proof of the 
above result could be carried out. In order to establish the equivalence between the 
action functional derived here and the action functional for the microscopic process 



derived in [33] we should think (3.37) as a maximization problem and use the calcu- 



lus of variation theory. Thus denoting H[x^t) the maximizer of (3.37), we have by 
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525 definition that for any appropriate test function <l> 



526 



543 



which can be written as 



- ^ (3.39) 



527 Substituting (3.39) into (3.37) follows that 

A(^,^,) =< d,H,L[^]d,H >L2 (3.40) 

528 and thus the rate function equals in the limit to 



Jo Jo 



5ot(^)=/ / L[md,Hydxdt (3.41) 



529 which is exactly the microscopic action functional given by (2.33). 

530 Remark: While for exchange (i.e. Kawasaki) dynamics the action functionals be- 

531 tween the Langevin approximations and the underlying microscopic process are asymp- 

532 totically equivalent, this is not true for adsorption/desorption (i.e. Glauber) dynam- 

533 ics. Indeed, both Langevin approximation [30 , [20 and Hanggi correction [34 result 

534 in action functional which are asymptotically different from the action functional of 

535 the underlying microscopic process derived in [HI p. 146]. Moreover, the action func- 

536 tionals of an SDE driven process is generally of weighted quadratic form [S^ while 

537 the action functional of the microscopic adsorption/desorption process is far more 

538 complex. However, using as a starting point for Langevin approximation the dis- 

539 cretization of the microscopic action functional -similar to what we did in this paper- 

540 there might be a way to construct accurate Langevin approximations whose action 

541 functionals are asymptotically equivalent to the microscopic adsorption/desorption 



542 process. 



4. Numerical Results. The objective of this Section is to study pattern forma- 

544 tion in surface diffusion using the proposed Langevin-type models. Since the stochastic 

545 fluctuations of the proposed models are directly derived from the microscopic process, 

546 the exploration of the pattern morphologies on the complex energy landscape of the 

547 particle system is well emerged. Moreover, authors consider that it is important to 

548 promote reproducible research hence the code written for the production of the flgures 

549 as well extended benchmark simulations is available online and it can be found at 

550 www .math .umass . edu/^pantazis/ source/patternForination_FigsCode . zip 



551 4.1. Numerical Schemes. In the previous Section space discretization (i.e. 

552 semi-discretization) was considered in detail. The flnal step in order to simulate the 

553 derived models on computers is to discretize the time, too. Since our primal goal 

554 is to highlight the space discretization, we keep the time discretization as simple 

555 as possible. Thus, a simple predictor-corrector (PC) Euler scheme which has 1st 

556 order weak convergence |23] is suggested. Of course implicit schemes or higher order 
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schemes such as Milstein's or derivative-free Runge-Kutta method could be used, 
however, they are computationaUy expensive especiaUy for high dimensional systems 
such as the studied. 

In order to highlight the implementation details, we restrict without loss of gen- 
erality only to DLM. Then the PC Euler scheme at n-th iteration is given in matrix 
form by 



Xn+1 =Xn^ [(1 - a)u{Xn) + au{Xn)]M + v{Xn)/^Wn 



(4.1) 



where At is the time-step while /SWn is a vector of independent zero-mean Gaussians 
with covariance matrix At/. Initial value of the lattice configuration denoted by Xq is 
also given while a is a weight factor which we set to 0.5 (trapezoidal rule). Since the 
size of the matrix v is x even though only + 1 of its diagonals are nonzero, 
it cannot be represented as a matrix in a computer memory hence we rewrite it - 
as well the drift term- in a compact implementable representation. For the general 
d-dimensional case, assume that k = {ki,...,kd) is a multi-index that denotes the 
position of the k-th. variable an is the unitary vector with 1 at position i. Then, 
the k-th. element of the drift term is given by 



-(^/c+ei(^n) -\- Lk{Xn)){Fk^ei{^n) - Fk^^n)) 



+ l^{Lk{Xn) -\- Lk-ei{Xn)){Fk{Xn) — F/c_e.(X^)) 



(4.2) 



572 where Fk{X) = —j3U{k^ X) + log ^^xlk) w^iile the k-th element of the stochastic term 

573 is given by 

VkAXn)MVn{l)=Y, 



d r 



1=1 



(4.3) 



^(Lfe(X„) + Lk-eM)AWi,{k - Bi) 



where ~ iV(0, At/^d) is a zero-mean Gaussian vector while and W^, are 
independent random vectors. 

In time discretization, similarly to space discretization, there are issues to be 
resolved. One such crucial issue is the choice of time step. At, which here were chosen 
heuristically using the following rule 



(4.4) 



keCrr 



which means that the average difference of the process in one step is controlled by 
S. After many experiments on a large parameter regime, we set S = 10~^ which 
is a compromise between stability and efficiency of the algorithm. Another artifact 
of time discretization is that the probability of X^+i leaving the admissible domain 
[0, 1]^ is 1 making the algorithm to diverge. A simple solution to this problem is 
that whenever there is a element of X^+i outside [0, 1] then the stochastic term is 
eliminated and only the drift term is considered. This is enough since the drift term 
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"push back" the value in the admissible interval. However, the cost to be paid is that 
we introduce bias which is proportional to the times the process leaves the admissible 
domain which of course depends on the time step, At. In our simulations, due to the 
specific choice of time step, the percentage of hitting the boundary values was less 
than 0.01%. 

4.1.1. Sources of CPU Acceleration. The most time-consuming part of the 
numerical algorithm is the computation of the potential U{k^Xn) at each step for 
all k G jCrn- This function is actually the convolution between the CG interaction 
potential and the lattice configuration. Thus an efficient method for computing the 
convolution between two function is through Fourier transform. Indeed, it holds that 



U{Xn) = J * = J'-'UiOMO} (4.5) 

597 where denotes the inverse Fourier transform while J(^) and the Fourier 

598 transforms of J and X^, respectively. 

599 Using multiplication in Fourier space instead of convolution in physical space 

600 makes the proposed method eventually independent of the interaction length. Indeed, 

601 the computational cost of one step of the numerical SDE solver is dropped from 
0{M^{L/qY) to 0(M^ log M^). Thus a huge computational gain is achieved for long 
range or mid range interaction potentials. This computational gain is a tremendous 

604 difference between the SDE approximations and the null event CGMC method which 

605 stems from the fact that in an SDE step the potential of all cells is needed while in a 

606 CGMC step the potential of only one cell is incorporated. Finally, the computation 

607 of convolution in Fourier space relates the proposed finite-difference method to the 

608 (pseudo-)spectral methods at least as concerns the computational cost. 
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4.2. Linear Stability Analysis. One fast and standard approach to roughly 
explore the behavior of the diffusive particle system at different parameter regimes is 
linear stability analysis of the mesoscopic PDF [36 . In connection with Fig. |2.1[ a), 
linearized techniques belong to the mean-field class of models where most of the 
atomistic details have been integrated out. Generally, linear stability analysis identi- 
fies when a spatial perturbation added to a uniform solution of the PDF would either 
eliminate of grow in time [37], [38]. Thus if we disturb a constant solution of the 
mesoscopic PDF 



aip = V-|L[p]V^| (4.6) 

by a spatially periodic perturbation of the form e^^e^^^ then the dispersion relation 
between the perturbation growth rate. A, and wavelength or mode, ^, is given by 



X^ = Tom\'L[co] 



co(l - Co) 



(4.7) 



where J(-) is the 2D Fourier transform of the continuous-space interaction potential, 

620 Co is the mean coverage and L[cq\ is the mobility either of Metropolis or Arrhenius 

621 dynamics for the constant density function p{t^x) = cq. 
In order to observe phase transition phenomena -in our case pattern formation- 



there should exist positive growth rates. From (4.7) we could predict that phase 
624 transitions occur when there exists at least one wavenumber ^' such that PJ{i') > 
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co(i-co) • Moreover, we could also predict from the same relation the most prominent 
size of the patterns. Indeed, the wavelength with the largest growth rate which is 
the wavelength that maximize the Fourier transform of the interaction potential (i.e., 
^max = argmax^ J{(,)) should dominate. Even though the following Section takes 
into account the information gained from linear stability analysis, it also reveals its 
limitations especially at critical parameter regimes. 

4.3. Pattern Formation Simulations. In order to perform reliable bench- 
mark simulations, it is necessary to utilize medium to large lattice domains. However, 
CGMC algorithm is prohibitively slow for large lattices resulting in the inability of pro- 
viding sufficient statistics for comparison. Thus we perform limited benchmark simu- 
lations and relied on the theoretical results obtained in previous Sections. Neverthe- 
less, we present the CPU time comparisons between CGMC algorithm and Langevin 



approximations. Table 4.1 shows the CPU execution time for null event CGMC algo- 
rithm and PLM2 model for Arrhenius dynamics. We prefer PLM2 model because it 
is the model with the most CPU-demanding (see (3.31)) among the Langevin mod- 
els. It is evident from the Table that PLM2 scales linearly as the size of the lattice 
is increased while CGMC scales super-linearly due to the fact that the time step in 
CGMC is inverse proportional to the lattice size. Moreover, PLM2 is about 10-20 
times faster from CGMC algorithm for relatively large lattices {N = 2^) achieving a 
significant time acceleration. 





CGMC 


PLM2 


7V = 2^ g = 2^ 


1.5 X 10^ 


3.3 X 10^ 


N = 2^,q = 2'^ 


3.6 X 10"^ 


2.5 X 10^ 



Table 4.1: CPU execution time in seconds of null event CGMC and PLM2 model for 
Arrhenius dynamics. Both algorithms run until final time T = 100. For N = 2^ both 
algorithms have converge to equilibrium while they have not for N = 2^. 



Proceeding now to the study of pattern formation phenomena in surface diffu- 
sion, an appropriate interaction potential should be chosen. Following [39] and [36], 
patterns are formed when interaction potential is attractive at short range resulting 
in microphase separation and repulsive at long range so as they do not coalescence. 
A typical choice of attractive/repulsive interaction potential is Morse potential given 
in a general form by 

T/ 7 \ ^1 \\^-y\\\ JiXi [ \\x-y\\ 

(4.8) 

where Ji is the potential strength while r^,! and r^,! are the attractive and repulsive 
length scales. Note that in order to have short range attractive and long range re- 
pulsive interaction potential it should hold r^,! > r^,!- The ratio between attractive 
and repulsive forces is determined by the repulsion strength, xi- The 2D Fourier 
transform of Morse potential is 



656 hence based on linear stability analysis the most prominent wavelength is the maxi- 
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657 mum of the Fourier transform of the interaction potential given by 



iiCTi=-^i/ y^^^ (4.10) 

658 where Ri = > 1 and it should hold 1 < y/xiRi < Ri so as a real- valued dominant 

659 mode is obtained. Moreover, the rate of growth of the dominant pattern size which 

660 is crucially determined from the value of the interaction potential at mode ^f^^^ (see 

661 (|4.7[)) equals 



i.«r") = ..''i^(i-,^) (4.11) 

However, the decay of the Fourier transform of the interaction potential is of poly- 
nomial order which is slow and under the presence of stochastic fluctuations patterns 
are irregular. In order to obtain nearly periodic configurations another interaction 
potential which is also called Morse potential should be utilized. In recent years, 
this potential had been applied for the study of pattern formation [40], [12] and it is 
defined as the difference of two Gaussian kernels, i.e. 

T f T\ -^2 f \\x-y\f\ J2X2 ( \\x-y\\^ \ 

u--y'^x2.wra^j^) '= ^-p [-^^ j " ^'"p ) 

'(4.12) 

where, similar to previous interaction potential, J2 is the potential strength, X2 is 
the repulsion strength while ra,2 and rr,2 are dimensionless length scales or attraction 
and repulsion, respectively. The 2D Fourier transform of this variant of the Morse 
potential is given by 

^2(0 = J2 exp I I - J2X2 exp I I (4.13) 

which is again a difference of two Gaussian kernels. Notice that the decay rate of 
the Fourier modes are now exponential. Since our primal interest is to produce con- 
figurations of patterns which are stable and nearly periodic we present most of our 
results using J2(-)- Moreover, the most prominent size of the patterns is related to 
the maximum value of the Fourier transform of J2(-) and it is obtained at 




\\^T''\\ = —\ hr. (4-14) 



677 where R2 = -r^ > 1 is the repulsive to attractive ratio while it should hold X2^2 > 1- 

678 The growth rate of the dominant wavelength is given by 

M^rn = J2(xi?2)^(l - (4.15) 

679 The study of pattern formation is performed using the variant of Morse potential. 



J2(-). Fig. 4.1 shows configurations of the system at equilibrium for various parameter 
values. Specifically, the size of the lattice is N = 2^ while the coarsening factor is 
q = 4. Interaction strength is J2 = 1 with inverse temperature is /3 = 12. Attraction 
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683 and repulsion length-scales are set to ra,2 = 5 and rr,2 = 10, respectively, while two 

684 different repulsion strengths, X2 = 0-4 (left column) and X2 = 0-8 (right column) are 
applied. Arrhenius dynamics with diffusion rate d(3 = 0.27 is used for this simulation 



685 



686 while the preferred numerical scheme was PLM2 with step size suitably chosen for 
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687 each case such that (4.4) is approximately valid. 

688 Based on linear stability analysis, we expect that patterns do occur for Fig |4.1[ a),(c), 

689 (d)&(e) but not for Fig 4TJb)&(f) because the growth rate as it is calculated from the 



690 dispersion relation (4.7) is negative for all modes. However as it is evident from the 

691 figure, patterns are formed in any case. Of course, patterns in Fig |4.1[ b)fc(f) are much 

692 more noisy exactly due to the fact that the growth rate of the dominant wavelength 

693 is (positive but) very small. Furthermore, the stochastic fluctuations of the model are 

694 important since patterns with different sizes are observed in each configuration. This 

695 result is in accordance with the CGMC runs performed in [36^ and it is far from the 

696 configurations obtained when deterministic models [12] were used where patterns are 

697 almost uniform. Additionally, changing the mean coverage, cq, dots, labyrinths and 

698 inverted dots are observed. Similar experimental images were shown in [27^ where 

699 surface diffusion of lead (Pt) on a copper (Cu) layer were studied. A final observation 

700 is that looking at the two columns of the Figure, the size of the patterns is decreased 



as repulsion strength, X2, is increased as it is expected from (4.14) since the dom- 
inant size of the patterns is inverse proportional to the wavelength. Intuitively, it 
can be also explained by the fact that strong long range repulsion leads to even less 

704 coalescence of patterns as time evolves. 

705 The final numerical experiment of this section is the comparison of the two at- 

706 tractive/repulsive interaction potentials Ji(-) and J2(-)- The motivation for this ex- 
periment stems from the fact that even though the prominent size of the patterns 
are chosen to be equal -based on linear stability analysis- for both potentials, the 
behavior of the overall system is expected to be different due to the different decay 
of the modes. As already stated, the decay of the modes for Ji(-) is polynomial while 
the decay of the modes for J2(-) is exponential hence we expect that the use of Ji(-) 
will produce a richer class of patterns making for instance the control of the size a 
rather difficult task. The configurations obtained at equilibrium using the Morse po- 

714 tential Ji(-) as well its variant J2(-) are shown in Fig. |4.2[ a) and (b), respectively. 
The control parameters of the interaction potentials was appropriately chosen so as 
the dominant modes of the Fourier transform be equal (i.e. ^^^^ = ^2^^^ = 0.15) as 
well their growth rates be equal (i.e. Ji(^i^^^) = ^2(^2^^^) = 1)- In order to specify 
all the parameters of the interaction potentials we further set Xi = X2 = 0.5 and 
Ri = R2 = 4: while the remaining parameters of the system are set to /3 = 10 and 
df3 = 0.5. 

721 By visual inspection of Fig. |4.2[ it can be stated that the distribution of the sizes of 

722 the patterns is more diverse for the original Morse potential compared to its variant. 

723 Moreover, as the histograms of the radius of the patterns suggest, the prominent 

724 radius in both potentials is 6 lattice sites which is comparable to the expected radius 

725 of the patterns being in this case = = 6.6 lattice sites. Finally, one simple 

726 approach to quantify the diversity of the pattern sizes is to compute the standard 

727 deviation of the radius of the patterns which is 2.4 and 2.0 for the original Morse 

728 potential and its variant, respectively. Overall, as it was predicted by linear stability 

729 analysis, original Morse potential produce a larger class of pattern sizes compared to 

730 its variant. 
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(a) CO = 0.2 & X2 = 0.4 



(b) CO = 0.2 & X2 = 0.8 
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(c) CO = 0.5 & X2 = 0.4 



(d) CO = 0.5 & X2 = 0.8 
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(e) CO = 0.8 & X2 = 0.4 



(f) CO = 0.8 & X2 = 0.8 



Fig. 4.1: A huge variety of patterns (dots, labyrinths, inverted dots) are produced 

at different parameter regimes corresponding to the complex landscape created by 

the competing interactions and the various conserved surface coverages, cq. Also, 

quantities such as the size of the patterns can be controlled by the system's parameters. 
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Configuration Radius Histogram 

(b) Applying J2 of \4.12\ 

Fig. 4.2: Using different attractive/repulsive interaction potentials configurations with 
different characteristics are obtained. 



731 5. Conclusions. In this paper, we derived models which served as approxima- 

732 tions of the CG process for the study of pattern formation on surfaces. Our starting 

733 point was an appropriate space discretization of the SPDE which lead to a system 

734 of SDE of Langevin type. Inspired by both the microscopic level and the mesoscopic 

735 level, the proposed models inherit properties from both levels. Indeed, (a) finite time 

736 estimates on the weak error between CG process and the process driven by the pro- 

737 posed models were obtained, (b) we showed that the action functionals between the 

738 microscopic model and the proposed are asymptotically equivalent which is a direct 

739 consequence of the fact that the proposed models are a direct discretization of the 

740 action functional (i.e. of the SPDE) and (c) by a perturbation of order O(^) either to 

741 the drift or the invariant measure, the derived models satisfied DBG hence the invari- 

742 ant measure of the approximation process is known. Hence the derived approximation 

743 models control the error at finite, long and infinite time scales. 

744 Additionally, the knowledge of the invariant measure revealed a very interesting 

745 observation -to our best knowledge never stated before- which says that the space 
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discretization of the SPDE for Metropolis dynamics lead to a system whose tempera- 
ture was perturbed by a factor of K This observation asserts that the discretization 
of a SPDE may produce artifacts and bias to the numerical results when q is small. 
Moreover, increasing or decreasing the power of the noise, which is straightforward 
for the suggested models by suitably scaling of the order parameters q and m, we are 
able to zoom in or out to more or less atomistic details of the system. In connec- 



tion with Fig. 2.1 'a), increasing or decreasing the power of the noise results in the 
translation of the models to the left towards microscopic level or to the right towards 
mesoscopic level, respectively. The controlled-error approximation and the microscop- 
ically derived fluctuations allow us to view the proposed models as "bridges" between 
molecular and continuum (S)PDE models of diffusion processes. Based on this reli- 
able intermediate models, it may be possible to consider hybrid micro/macro models 
bridging the gap between algorithms with different spatial scales. We also refer to 
recent work in related hybrid models in fluctuation hydrodynamics [41 . 

Finally, as concerns the study of pattern formation phenomena through a self- 
assembly mechanism, we efficiently reproduce the sizes and types of patterns ex- 
perimentally observed in previous studies [23 . The role of noise is critical for the 
systematic exploration of the complex energy landscape of the system. As it was 
evident from the Figures, the choice of the interaction potential as well the variation 
of the system's parameters significantly affects the size and the shape of the patterns. 
Additionally, having the invariant measure of the process, one of our next goals is to 
perform sensitivity analysis using the method developed by Majda and Gershgorin 
[22] which exploits the Fisher information at equilibrium. Furthermore, another im- 
portant application we are interested in is the control of the pattern's properties. By 
varying the parameters of the system such as the mean concentration or the temper- 
ature or the repulsion strength in a controlled way we will be able to design patterns 
with specified shapes, sizes or even orientations. However, in order to perform op- 
timal control we need appropriate mesoscopic observables for the patterns which is 
also under our research investigation. Particularly, defining appropriate mesoscopic 
observables using tools from image processing (see right column of Fig. 4.2 as a pre- 



liminary example of such tools) and pattern recognition is one of our immediate goals. 
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Appe ndix A. Asvmptotics . In order to compute the forma l asymptotics of the 
drift term ( |3.2|) , diffusion matrix (3.5) and correction term (3.24) for both dynamics 
we need theToTlowing Taylor series expansions. IVIobility is expanding up to the second 
order derivative given by 

Lk±i(p) = Lk(p) ± —dxLk(p) + -^dxxLk(p) + O(^) 

while potential is similarly expanded up to third order derivative by 

U{k±l,p) = U{k,p)±-dxU{k,p) + -^dxxU{k,p)±-^dxxxU 

Finally, the difference of the logarithms which stems from the entropy term is ex- 
panded as 
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and similarly for log(l - pk) - log(l - pk±i)' 

Now, we are able to compute the formal asymptotic for the drift term as 

uk{p) = liLk+lip) + Lk(p)){-(3[U(k + l,p)- U(k,p)] + [log(p(xfc+i)) - \og(p(xk))] - [log(l - p(xk+i)) - log(l - p(xk))]} 
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870 Similarly, the weak form of the diffusion matrix (i.e. covariance matrix) has the 

871 following formal asymptotic 



\ k,j l,i / k,l 

= ^ [Dk,k-iH{xk-i) + Dk,k(p2{xk) + Dk,k+iH{xk+i)] (t>i{xk) 

k 

= X] [Dk,k-l{(f>2{xk-l) - (f)2(xk))((f)l(xk-l) - (f)l{xk))] 
k 

= ~yZ K^k-l + Lk)((t)2(xk-l) - (t)2(xk))((t)l(xk-l) - (t>l(xk))] 
^k 

= -^Lk [((f)2(xk-l) - (f)2(xk))((t>lixk-l) - (t>lixk)) + ((f)2(xk + l) - (f)2(xk))i(l>lixk + l) - (/>l(xfc))] 



2 
qm 



-dx(t)i{xk)dxct)2{xk) + 0(l/m'^) 



j L(p(x))d:,4>i(x)d:,4>2(x)dx + 0(l/m^) 



(A.2) 



873 Finally, the asymptotic expansion of the "correction" term for Arrhenius dynamics 

874 can be alternatively written as 



Ckip) = ^^xx {(1 - 2p{xk))eM-f3U{k,p)) - i3q(J(0) + J(l))Lfc(p)} + O(^) 

= {-2d:,p(xk) exp(-/3t7(fc, p)) - 13(1 - 2p(xk))d:,U(k, p) exp(-/5[7(fc, p)) 

- /3g(J(0) + J(l))[(l - 2p{xk))d:,p{xk) exp(-/3f7(fc, p)) - (3da,U{k, p)Lfe(p)]} + O(^) 
2dxp(xk) , ^2 



rn^ I Vp(a:fc)(l - P{xk)) 
/5— ^^^^^^%^(g(J(0) + J(l))a.p(xfc) + d.U{Kp))) Lk{p)\ + O(^) 



(A.3) 



876 Appendix B. Detailed Balance Condition and Disc rete Free Energy 

877 Decrease. The computation of the inner products shown in (3.16) is given next. 

878 After an integration by parts and taking advantage of the periodic boundary condition 

30 



879 we obtain that 

Z <Mf,g>L2^^^ 



2? J dpk+i \dpk J dpk-i \dpk 



=--f 

2qJ 



k+1 



Opk 



dpk+1 



J2 {Lk+Lk-i)-^e 



ke£n 



dpk 



dpk 



2q J 

E (A 



9/ .-,B(p)i?SL 



dpk 



2qJ ^f-f dpk 



99 , , df 



kec 

-qE{p) 



(B.l) 



881 Discrete free energy functional, E{p)^ is decreasing over time. Indeed, taking once 

882 again advantage of the periodic boundary condition we obtain 



dt 



dpk dt 



1 ^ dEjp) 



{Lk+i + Lk) 



dE{p) dE{p) 



keCr, 



dpk+i dpk 
dE(p) dE{p) 



■ {Lk + Lk-i) 



dE{p) dE{p) 



dpk dpk-i 



dpk+1 dpk 
dE{p) dE{p) 



dEjp) 

dpk 
dEjp) 

dpk [dpk+i dpk J 2 



i J2 (Lk+Lk-i 



keCr, 



dEjp) 

dpk 
dEjp) 
dpk+i 



dE{p) dE{p) 



dpk dpk-i] 
dE{p) dE{p) 



dpk 



+1 



= E (^fe + l + ^k) 



dE{p) dE{p) 



dpk^ 



dpk . 



< 



dpk 



(B.2) 



883 since mobihty is always a non-negative function. 
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